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Abstract 

An automated flutter boundary tracking 
procedure Is presented for the efficient 
calculation of transonic flutter boundaries. 
The new procedure uses aeroelastlc responses to 
march along the boundary by taking steps In 
speed and Mach number, thereby reducing the 
number of response calculations previously 
required to determine a transonic flutter 
boundary. The tracking procedure reduces 
computational costs since only two response 
calculations are required per Mach number and 
provides a complete boundary in a single job 
submission. Flutter boundary results are 
presented for a typical airfoil section 
oscillating with pitch and plunge degrees of 
freedom. These transonic flutter boundaries are 
In good agreement with "exact" boundaries 
calculated using the conventional time-marching 
method. The tracking procedure was also 
extended to include static aeroelastlc twist as 
a simulation of the static deformation of a wing 
and thus contains all of the essential features 
that are required to apply It to practical 
three-dimensional cases. Application of the 
procedure Is also made to flutter boundaries as 
a function of structural parameters, the 
capability of which Is useful as a design tool. 


Nomenclature 

a nondimenslonal elastic axis location, 

positive aft of midchord 
b airfoil semichord 

c airfoil chord 

c* lift coefficient 

c» steady-state lift coefficient 

c m ° moment coefficient about pitching axis 

c m steady-state moment coefficient 

about pitching axis 

h plunge displacement, positive down 

[K] structural stiffness matrix 

K a pitch spring constant 

m airfoil mass per unit span 

M freestream Mach number 

[M] structural mass matrix 

r a airfoil radius of gyration about 

elastic axis 


r C g airfoil radius of gyration about 

center of gravity 

s Laplace transform variable, o+lu 

t time In seconds 

U freestream velocity 

Uf flutter speed 

V speed index, U/(bui a /IT) 

Vf flutter speed index, Uf/(bw a /H) 

x distance aft of leading edge 

x a nondimenslonal distance from elastic 

axis to mass center 

a airfoil angle of attack, positive 

leading edge up 

oq airfoil mean angle of attack 

a r wing root angle of attack 

5 dominant damping of aeroelastlc 

response, positive for stable response 
u airfoil mass ratio, m/npb? 

p freestream air density 

a damping coefficient of Laplace transform 

variable 

oi frequency of aeroelastlc response 

uf, uncoupled natural frequency of 

plunge 

oi a uncoupled natural frequency of 

pitch about the elastic axis 


Introduction 

Flutter is frequently a limiting factor in 
the performance of aircraft In transonic 
flight. Because of this limitation, It Is 
highly desirable to be able to accurately 
predict transonic flutter characteristics, 
especially during aircraft design. Furthermore, 
since design Is an Iterative process, it Is 
Important to be able to perform transonic 
flutter calculations inexpensively. 

Research on developing nonlinear transonic 
flutter prediction techniques has progressed 
with the development of transonic small-dls- 
turbance (TSD) computer codes. The nonlinear 
techniques typically require the calculation of 
time responses to determine the stability of the 
aeroelastlc system. In the time response 
analysis, the structural equations of motion are 
coupled to transonic aerodynamic codes using a 
numerical Integration procedure. For example, 
Ballhaus and Goorjian 1 first reported the cal- 
culation of transonic aeroelastlc responses for 
an oscillating airfoil with a single pitching 
degree-of -freedom (dof). Time responses were 
presented for the NACA 64A006 airfoil at M = 
0.88. The calculations were performed by 
simultaneously Integrating the structural 
equation of motion with the unsteady aerodynamic 
solution procedure of their LTRAN2 2 code. 
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The time-marching aeroelastic solution 
technique has subsequently been refined and 
extended by a number of researchers. Guruswamy 
and Yang 3 performed aeroelastic response 
calculations for a plunging and pitching NACA 
64A006 airfoil at M » 0.85. Flutter speeds 
selected from a separate flutter analysis Indeed 
resulted In neutrally stable time responses. 
Time-marching aeroelastic responses for airfoils 
Including an aileron pitching dof were presented 
In Refs. 4-6. Rlzzetta 4 demonstrated sig- 
nificant nonlinear effects using relatively 
large Initial conditions. Yang and Chen 5 and 
Yang and Batina 6 presented three dof time- 
response results for the NACA 64A006, NACA 
64A010, and MBB-A3 airfoils. Edwards, et al. 7 
presented nonlinear transonic flutter boundaries 
Including angle-of-attack effects. The results 
were obtained using a more accurate structural 
Integrator which was based on the state- 
transition matrix. More recently, Batina and 
Yang 8 have Investigated the effects of active 
controls on transonic aeroelastic responses, and 
subsequently assessed the accuracy of state- 
space aeroelastic modeling. Berry, Batina, and 
Yang 9 have Included viscous effects In the 
calculation of transonic aeroelastic responses 
by using an Integral boundary layer routine 
coupled to an Invlscld TSD code. Furthermore, 
Borland and Rlzzetta 10 have applied the time- 
marching aeroelastic response technique to a 
three-dimensional case. In Ref. 10, nonlinear 
transonic flutter results were reported for a 
rectangular wing with parabolic arc airfoil 
section. 

Although the time-marching aeroelastic 
procedures are relatively well developed for the 
determination of transonic flutter boundaries, 
the method can require large amounts of computer 
resources. In general, several aeroelastic 
responses need to be calculated for a range of 
flight speeds at a given Mach number. The 
flutter point is then estimated by Interpolation 
of the damping from stable and unstable 
responses. To determine the flutter boundary, 
the procedure Is repeated for each Mach number 
of Interest. The repetitive nature of the 
time-marching aeroelastic method (hereafter 
referred to as the conventional method), can be 
computationally quite expensive as well as 
manpower Intensive. The calculation of a 
transonic flutter boundary, for example, can 
require as much as two to four weeks of elapsed 
analysis time. 

What Is needed Is a second generation 
method, which Is less expensive and time 
consuming. Therefore, It Is the purpose of this 
paper to present the development of a transonic 
flutter boundary tracking procedure, which 
significantly reduces the elapsed analysis time 
and lowers computational costs. The new 
procedure uses a time-marching aeroelastic 
response code to systematically march along the 
boundary by taking steps In speed and Mach 
number. This tracking reduces the number of 
aeroelastic response calculations which are 
required for the conventional method. The 
objectives of this research were: (1) to 
develop an efficient automated flutter boundary 
tracking procedure which utilizes a time- 
marching aeroelastic response code, (2) to 
verify ’ the procedure by making detailed 


comparisons with conventional flutter solutions, 
(3) to Improve the procedure by including static 
aeroelastic twist, and (4) to demonstrate the 
robustness and utility of the new procedure by 
applying It to a variety of aeroelastic cases. 

In this study, the tracking procedure Is 
applied to a simple aeroelastic system 
consisting of a typical airfoil section 
oscillating with pitch and plunge degrees of 
freedom. The procedure Is coupled to the time- 
marching aeroelastic response analysis within 
the XTRAN2L 11 unsteady transonic small - 
disturbance code. The paper presents a detailed 
description of the flutter boundary tracking 
procedure along with representative results and 
comparisons which assess the new method. 


Computational Procedures 

In this section, the time-marching 
aeroelastic method of calculating transonic 
flutter speeds Is first briefly described. A 
more detailed description Is given of the 
flutter boundary tracking procedure Including 
the step-by-step numerical Implementation. 
Finally, the Inclusion of static aeroelastic 
twist Into the calculations to simulate static 
deformation under load Is also described. 

Time-Marching Aeroelastic Solutions 

The equations of motion for a typical 
airfoil section with pitch and plunge degrees of 
freedom may be written In matrix form as 7 


CeT 


[M]. 


} * [K] 




0,1 
r j 


I ,u,‘ 

Tip 'F' 


2(C-C m >( 
m m o 'J 


(1) 


where [M] Is the mass matrix, [K] Is the 
stiffness matrix. and the dot denotes 
differentiation with respect to time t. Plunge 
and pitch displacements h and a, respectively, 
are perturbations about the mean airfoil 
position. The lift and moment coefficients 
c 4 and % , respectively, are the 
steady-state values calculated at the mean angle 
of attack oq , which are subtracted from the 
total coefficients. In the time-marching 
solution, the equations of motion are coupled to 
an unsteady aerodynamic code for 
time-integration. In this study, Eq. (1) Is 
simultaneously Integrated with the aerodynamic 
solution procedure of the XTRAN2L code using the 
modified state-transition matrix Integrator of 
Edwards, et al. 7 Details of the solution 
procedure may be found In Ref. 7. An Initial 
plunge displacement Is typically used to start 
the Integration of Eq. (1) by setting h(0) = 
0.01, and a(0) = fi(0) = a(0) = 0. 


In the conventional method, the flutter 
analyst must repetitively calculate time- 
marching aeroelastic responses for a range of 
speeds U, at a given Mach number. Generally 
three to four speeds are selected for response 
calculations to determine one stable and one 
unstable flight speed. The flutter speed Uf, 
Is calculated by Interpolating the dominant 
damping of the responses to the speed where 
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damping is zero. The procedure is then repeated 
for each Mach number of Interest. 

Damping of the aeroelastic responses Is 
estimated from the response curves using the 
modal identification technique of Bennett and 
Desmarais. 12 The modal estimates are determined 
by a least squares curve-fit of the 
responses using complex exponential functions of 
the form 


X(t) - A + I e° / 
J=1 


[AjCOsfu t) + BjSin^t)] (2) 


where m is the number of modes selected here to 
be equal to two. The damping of each mode Is 
calculated using 


c 


j 



(3) 


The dominant damping of the aeroelastic system 
is determined by selecting the smallest value of 
Ci. As an example, the top part of Fig. 1 
snows a typical pitch response and modal 
curve-fit. In the lower part of Fig. 1 are the 
two component modes determined from the modal 
fit. As shown in Fig. 1, the lower frequency, 
constant amplitude mode is dominant, and 
therefore determines the stability of the 
aeroelastic system. 


make decisions on continuing the solution 
procedure. Therefore, the method can be quite 
manpower Intensive with the calculation of a 
transonic flutter boundary requiring as much as 
two to four weeks of elapsed analysis time. 

Flutter Boundary Tracking Procedure 

The flutter boundary tracking procedure is 
an algorithm for calculating a transonic flutter 
boundary in an efficient and automated fashion. 
In this study, the procedure is coupled to the 
time-marching aeroelastic response analysis 
within the XTRAN2L code. However, the procedure 
is general enough to be coupled to any 
aeroelastic response code. 

Starting conditions are required to begin 
the flutter boundary tracking procedure at Mach 
number Mj, as shown in Fig. 2(a). The 
starting conditions are obtained within the 
tracking procedure code using the conventional 
method. Several responses are computed for 
varying speed indices to determine values Vj 
and V? which bracket the flutter boundary. 
The first flutter speed index Vf. is 
calculated by interpolation of the dominant 
damping of these responses. The starting 
conditions for the tracking procedure at Mach 
number M} include the speed indicies Vj and 
V 2 , the dominant damping values of the 
responses at these speeds q and q>> 
respectively, and the flutter speed index 

Vf r 


In the conventional method, the modal 
identification code of Ref. 12 is normally run 
as a separate job step. Based upon the results 
of this program, the flutter analyst must then 
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Fig. 1 Modal curve-fit of aeroelastic pitch 
response and resulting component modes 
at flutter for the NACA 64A010A airfoil 
at M = 0.80. 


The procedure tracks the boundary from a 
given Mach number Mt to M 1+ i as shown in 
Fig. 2(b). The flutter speed index Vf is 

assumed to be known. The flutter speed index at 
M^+i is predicted by approximating the flutter 
boundary curve with the Taylor senes 



dV 


+ -i( M i + i- M i> 
dM 2 


+ I! B (M -1 - V 2 + 


(4) 


Since relatively small steps in Mach number are 
used, all terms of second or higher order are 
neglected. Since V = f(M) Is an implicit 
function specified by c(V,M) = 0, the slope of 
the flutter boundary can be expressed as dV/dM = 
-(3C/3M)/(3 c/ 3V). Rewriting Eq. (4) in delta 
notation gives 

al./ am 

V, • V ( "<« ‘ 151 


The change in damping terms are defined as 



AC V = ctVj.H,) - e(v 2 ,M i ) 

(6a) 

ac m - 

c( v 3 ,m, + i) - rtVj^), If V 3 = Vj 

(6b) 

ic M = 

c(v 3 , m i + i) - c(v 2 ,m 1 ), if V 3 = V 2 

(6c) 
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(a) 


V, 

(b) 

v i 

(c) 

Fig. 2 


O Response calculation 
O Interpolated flutter point 



calculation of first flutter speed 
Index Vfj and starting conditions. 


O Response calculation 
□ Predicted flutter point 
O Interpolated flutter point 



Mach number 


flutter boundary tracking from Mach 
number M< to Mf+i. 



continuation of tracking procedure 
over a range of Mach number values. 


Graphical description of flutter 
boundary tracking procedure. 


and the terms AV and AM are expressed as 


AV = Vj - V 2 

(7a) 

AM = M 1+1 - H, 

(7b) 


where c Is the damping of the dominant mode of 
the aeroelastlc response calculated for the 
given speed Index and Mach number combination. 
The speed Indices V^, V 2 , and V3 are the 
values at which response calculations are 
performed. Equations (5), (6), and (7) form the 
basis of the flutter boundary tracking 
procedure. 

Computationally, the procedure tracks the 
boundary from Mach number M$ to M 1+ j as 
described In steps 1 through 4. 

1. A converged steady flowfleld is calculated 
at conditions (Vf 1 , Mi ) to be used as 

the initial flowfleld for aeroelastlc 
response calculations at Mf+i. For the 
values of AM Investigated, tne numerical 
transient Induced by using an Initial 
flowfield at a neighboring Mach number is 
negligible. 

2. To determine the flutter speed index at 

Mach number Mi + i, aeroelastlc responses 
are calculated at the speed Index V3 
shown In Fig. 2(b). The dominant damping 
of these responses Is then determined for 
the evaluation of the Acm/ 4M term in Eq. 
(5). The Acy/AV term in Eq. (5) is 

calculated from Information already known 
at Mi. For the first prediction of 

Vf 1+1 from Vf.j with Eq. (5), the 

value of the speed index V3 is set equal 

to Vj. For all other predictions, the 
speed Index V3 Is set equal to Vj If 
the absolute value of (Vf^-V^) is 

less than the absolute value of 

( Vf -V 2 ). Otherwise, V3 is set 

equal to V 2 . The predicted flutter point 
determined using Eq. (5) Is labeled 
(Vf 1+1 )p In Fig. 2(b). 

3. Aeroelastlc response histories are 

calculated at V4 which is set equal 

to (Vf 1+1 ) p to determine the 

stability of the airfoil at the predicted 
flutter speed index. The dominant damping 
values obtained from the responses 
calculated at V3 and V4 are then used 

to calculate the final flutter speed index 
Vf i+1 by Interpolation or extrapolation 

of the damping. Figure 2(b) shows a 
typical interpolation case. 

4. The procedure is continued to the next Mach 

number by returning to step 1 for the 

calculation of a new converged steady 
flowfleld. When returning to step 1 the 
values for V lt V 2 , ci , C2. Vf , 

and Mi are set equal to V3, V4, 43, 

C4, Vf i+1 , and M 1+1 , respectively. 

By repeating steps 1 through 4 for a range 
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of Mach number values, a flutter boundary 
can be calculated as shown graphically in 
Fig. 2(c). 

The flutter speed index Vf 1+1 at Mach 

number M-j + i is calculated using linear 
interpolation or extrapolation of the dominant 
damping at speed indicies V 3 and V 4 in step 
3 of the tracking procedure. To demonstrate the 
validity of this assumption of linearity in 
damping for small changes in speed about the 
flutter point, the dominant damping values from 
several representative response calculations are 
plotted versus speed index in Fig. 3. The 
response histories were calculated for the NACA 
64A010A airfoil at M = 0.80, in Increments of 
speed index AV equal to 7% of the resulting 
flutter speed index. As shown in Fig. 3, the 
assumption of linearity is valid since the 
damping values (circles) lie very close to a 
straight line. 



Additionally, the flutter boundary tracking 
procedure is not constrained to determine 
flutter boundaries as a function of Mach 
number. The current implementation of the 
procedure allows for the calculation of flutter 
boundaries with respect to an aeroelastic 
parameter such as w, a, x a , r a , <^,, or 
u> a instead of Mach number. In this case, the 
mathematical formulation of the tracking 
procedure is the same as that described above 
with the chosen parameter being substituted for 
Mach number M. Computationally, the parameter 
and the corresponding step are chosen as input. 

Static Aeroelastic Twist 

To be able to conceptually apply the 
tracking procedure to realistic 

three-dimensional cases, the procedure was 
extended to include static aeroelastic twist as 
a two-dimensional simulation of the static 
deformation of a wing. The significance of 
static aeroelastic twist on typical section 
flutter at transonic speeds was demonstrated by 
Edwards, et al . 7 The static twist is determined 
by using a simple model of a linear root pitch 
spring and equating the aerodynamic pitching 
moment to the pitch spring restoring moment. 
The resulting static equilibrium equation may be 
written as 

K a<V“r> * > 2 (2b) 2c m (“ 0 .M) (8) 

0 


where (ag-a r ) is the static aeroelastic 
twist. In terms of a strip theory analysis, 
a r is the wing root angle of attack and ag 
is the local section mean angle of attack. The 
steady -state moment coefficient about the 
pitching axis c„, Is the same quantity that 
appears in Eq. °1). By nondimenslonalizlng, 
Eq. ( 8 ) may be rewritten as 


“o - a r + 5? C m<V M) 

* r a 0 


(9) 


which Is solved iteratively for the mean angle 
of attack ag. 


Fig. 3 Dominant damping vs. speed index for 
NACA 64A010A airfoil at M = 0.80. 


The mean angle of attack must be calculated 
to Include static aeroelastic effects in flutter 
boundary calculations. By modifying XTRAN2L to 
Include Eq. (9), a steady flowfield and mean 
angle of attack are obtained simultaneously. 
Specifically, the wing root angle of attack 
is fixed and ag is varied until a converged 
steady flowfield and mean angle of attack og 
are determined. For application to the con- 
ventional method, a static aeroelastic 
calculation should be performed before each 
response calculation. This adds the com- 
putational cost of performing one steady 
flowfield calculation per response to the total 
cost of a conventional method flutter solution. 
When including static aeroelastic twist in the 
flutter boundary tracking procedure, the mean 
angle of attack ag is determined during the 
steady flowfield calculation in step 1. For the 
cases presented a converged steady flowfield and 
mean angle of attack og are obtained in the 
same number of time steps required for the 
steady flowfield calculation without static 
twist. Consequently static twist is included in 
the tracking procedure at no additional cost. 
The mean angle of attack used for all response 
calculations at M 1+ i is set equal to the mean 
angle of attack calculated for flutter at M^ . 
For the values of AV and aM investigated, the 
changes in mean angle of attack ag when 
stepping from Mt to M 1+ j were small. 
Consequently, the static aeroelastic twist is 
lagged one Mach number for computational 
convenience. 


Results and Discussion 


Flutter boundary calculations were 
performed for a typical airfoil section 
oscillating with pitch and plunge degrees of 
freedom. The airfoils chosen were the NACA 
64A010 NASA Ames model (herein referred to as 
NACA 64A010A) and the MBB-A3 supercritical 
airfoil. The airfoil coordinates were taken 
from Ref. 13. Flutter boundary results are 
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presented for Case A of Isogai 14 * 15 which has 
normal modes similar to those of a streamwise 
section near the tip of a sweptback wing. The 
wind-off coupled plunge and pitch frequencies 
are 71.33 and 535.65 rad/sec, respectively. The 
pivotal point for the plunge mode is located 
1.44 chordlengths ahead of the leading edge. 
The pivotal point for the pitch mode is 0.068 
chordlengths forward of midchord. Specifically, 
the aeroelastic parameter values are a = -2.0, 
x a = 1.8, r a = 1.865, u = 60.0, = 100 
rad/sec, and io a = 100 rad/sec. 

Flutter boundary tracking results are 
compared with “exact" flutter boundaries 
calculated using the conventional method. These 
comparisons serve to verify the flutter boundary 
tracking concept as well as to assess the 
accuracy of the procedure. Flutter boundary 
tracking results are presented both with and 
without static aeroelastic twist. Static 
aeroelastic twist is included to further 
demonstrate the versatility of the procedure and 
to show the importance of static deformation on 
flutter behavior. Finally, flutter boundary 
tracking results are presented as a function of 
elastic axis location to show a further 
extension of the procedure. 

Flutter Boundary Tracking Results 

Calculations were performed for the NACA 
64A010A and MBB-A3 airfoils for Mach numbers 
ranging from M = 0.65 to M = 0.80. The mean 
angles of attack were a 0 = 1.0° for the NACA 
64A010A airfoil and a 0 = -0.5° for the MBB-A3 
airfoil. These angles of attack were chosen by 
Bland and Edwards 1 , for these airfoils, because 
they produce steady-state lift and shock 
locations that are approximately equal at the 
same Mach number. With lift and shock locations 
approximately equal, comparisons of flutter 
behavior can be made between the two airfoils. 
Steady pressure distributions for the NACA 
64A010A and MB8-A3 airfoils are shown in Figs. 
4(a) and 4(b), respectively, for Mach numbers 
from M = 0.75 to M = 0.80 in increments of 
0.01. For the NACA 64A010A airfoil (Fig. 4(a)), 
the flow is subcritical up to approximately M = 
0.76. Between M = 0.76 and 0.77, a shock wave 
forms on the upper surface near 45% chord. With 
increasing Mach number, the shock grows in 
strength and moves further downstream along the 
airfoil. The shock becomes relatively strong 
and is located near 68% chord at M = 0.80. For 
the MBB-A3 airfoil (Fig. 4(b)), the steady 
pressure distributions behave in much the same 
manner as for the NACA 64A010A airfoil. The 
upper surface shock wave forms between M = 0.77 
and 0.78, and grows in strength with increasing 
Mach number. At M = 0.80, the shock is located 
at 68% chord which is the same as that for the 
NACA 64A010A airfoil. The strength of the 
MBB-A3 shock wave, though, is considerably 
less. Furthermore, the pitching moment (about 
quarter chord) for the MBB-A3 airfoil Is much 
different in comparison with the NACA 64A010A 
airfoil because of the aft-loading. 

Transonic flutter boundaries for the NACA 
64A010A and MBB-A3 airfoils are shown in Figs. 
5(a) and 5(b), respectively. Flutter boundary 
tracking results are presented for steps in Mach 
number of = 0.01, 0.03, and 0.05. The three 



(a) NACA 64A010A airfoil at «o - 1-0°. 



Fig. 4 Steady pressure distributions for Mach 
numbers ranging from M = 0.75 to M = 
0.80. 


sets of results were obtained to assess the 
accuracy and robustness of the tracking 
procedure by making direct comparisons with 
conventional solutions. For the two airfoils, 
the flutter boundaries are quite similar. As 
pointed out in Ref. 16, the boundary for the 
MBB-A3 airfoil is nearly identical to that of 
the NACA 64A010A airfoil when it is shifted to 
the left by 0.01 Mach number. These flutter 
boundaries show the so-called transonic dip, but 
for the Mach number range considered, the 
boundaries do not define the minimum flutter 
speed. The Mach number range is restricted 
since potential codes are not reliable beyond 
about M = 0.80 for these airfoils unless the 
entropy generated by the shock waves is 

accounted for. 17 For the NACA 64A010A airfoil 
(Fig. 5(a)), all three sets of tracking results 
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Fig. 5 Comparison of flutter boundaries 

calculated with conventional and flutter 
boundary tracking procedures. 


are In very good agreement with the conventional 
solution. The average errors for AM = 0.01, 
0.03, and 0.05 are 0.39%, 0.74%, and 2.21%, 
respectively. The damping and speed index 
values that were used to determine the aM = 0.03 
tracking results are listed in Table 1. The 
conventional method flutter speed indicles are 

labeled Vf gxact in Table 1. The tabulated 

speed indicies indicate that the tracking 
procedure stays close to the boundary when 
marching from one Mach number to the next. For 
the MBB-A3 airfoil (Fig. 5(b)), the three sets 
of tracking results are in very good agreement 
with the conventional solution except at M = 
0.80. Here, the flutter points for AM = 0.03 
and 0.05 are slightly overpredicted because of 
the steepness of the boundary. The average 
errors for AM = 0.01, 0.03, and 0.05 are 0.46%, 
1.65%, and 4.34%, respectively. For the latter 
two cases, the average error is dominated by the 
overprediction of the flutter points at M = 
0.80. An option exists within the current 
Implementation of the flutter boundary tracking 
procedure to perform ^an additional response 
calculation per Macir number to improve 
accuracy. When this option was exercised for 
the MBB-A3 airfoil, the average errors for am = 
0.03 and 0.05 were significantly reduced to 
0.62% and 0.89%, respectively. The results 
presented demonstrate that the flutter boundary 
tracking procedure is accurate and robust. The 
tracking results were obtained with fewer 
response calculations than the conventional 
method results. A computational savings of 
approximately a factor of two was attained. In 
all of the remaining results to be presented, 
the step in Mach number was set equal to 0.02, 
and only two response calculations per Mach 
number were performed. 

Flutter Boundary Tracking Results Including 
Static Twist 

Flutter boundary calculations including 
static aeroelastic twist were performed for the 
NACA 64A010A and MBB-A3 airfoils. The flutter 
boundaries were calculated for Mach numbers 
ranging from M = 0.70 to M = 0.80. The wing 
root angles of attack were selected as a r = 
1.0° for the NACA 64A010A airfoil and a r = 
-0.5° for the MBB-A3 airfoil. These root angles 
of attack are identical to the mean angles of 
attack oq used in the previous section to 


Table 1. Flutter boundary tracking results for the NACA 64A010A airfoil with AM = 0.03. 


M 1 

V 1 

V 2 

c(Vj) 
X 10 3 

Eg 

M i + 1 

V 3 

■ 


c(v 4 ) 

X 10 3 


v f 

exact 






0.65 

1.64 

1.58 

- 5.95 

2.08 

1.59 

1.59 

0.65 

1.64 

1.58 

- 5.95 

2.08 

0.68 

1.64 

1.49 

-18.37 

0.99 

1.50 

1.50 

0.68 

1.64 

1.49 

-18.37 

0.99 

0.71 

1.49 

1.40 

-12.37 

-1.50 

1.39 

1.38 

0.71 

1.49 

1.40 

-12.37 

-1.50 

0.74 

1.40 

1.25 

-19.28 

-1.24 

1.24 

1.23 

0.74 

1.40 

1.25 

-19.28 

-1.24 

0.77 

1.25 

1.02 

-30.27 

-3.62 

0.99 

0.97 

0.77 

1.25 

1.02 

-30.27 

-3.62 

0.80 

1.02 

0.52 

-50.64 

5.37 

0.55 

0.55 
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allow for direct comparison of results with and 
without static twist. Steady pressure 
distributions for these airfoils at M = 0.80 are 
shown In Figs. 6(a) and 6(b). The speed 
Indlcles used for the steady pressure 
distributions computed with static twist were 
the resulting flutter speed values, Vf = 0.67 
and Vf = 0.94, for the NACA 64A010A and MBB-A3 
airfoils, respectively. For this case the 
static twist decreases the strength of the upper 
surface shocks and reduces the steady-state 
loading. For the NACA 64A010A airfoil (Fig. 
6(a)), the static twist lowers the mean angle of 
attack from = 1.0° to = 0.46°. This 
change In mean angle of attack weakens the shock 
and displaces It forward from 68% chord to 61% 


Without static twist 

With static twist 



(a) NACA 64A010A airfoil with otp - 1.0°. 


Without static twist 

With static twist 



(b) MBB-A3 airfoil with otp = -0.5°. 


Fig. 6 Steady pressure distributions at M = 
0.80, Including the effect of static 
aeroelastlc twist at flutter. 


chord. Additionally, a relatively weak shock 
forms on the lower surface at 49% chord. For 
the MBB-A3 airfoil (Fig. 6(b)), the static twist 
lowers the mean angle of attack from oq = 
-0.5° to oq = -1.5°. The affect of static 
twist on the pressure distributions of the 
MBB-A3 airfoil Is similar to that on the NACA 
64A010A airfoil. The upper surface shock wave 
Is significantly weakened and is displaced 
forward from 68% to 61% chord, which is the same 
shift In shock location that was determined for 
the NACA 64A010A airfoil. In general, static 
twist affected the steady pressure distributions 
of the two airfoils in a similar manner over the 
range of Mach numbers considered. The static 
twist angles for both airfoils are plotted as 
functions of Mach number in Fig. 7. These 
angles were computed at the speed indicies 
corresponding to flutter (to be presented 
subsequently). Results are only plotted from 
the flutter boundary tracking procedure since 
the conventional method twist angles are 
identical to plotting accuracy. As shown in 
Fig. 7, the twist angles decrease with 
increasing Mach number since the flutter speed 
index decreases. The static twist for the 
MBB-A3 airfoil is much larger (negatively) than 
that of the NACA 64A010A airfoil throughout the 
entire Mach number range considered. This Is 
attributed to the different pitching moment 
characteristics of the MBB-A3 airfoil due to the 
aft-loading. 

Flutter boundary tracking results, computed 
both with and without static twist, are shown in 
Figs. 8(a) and 8(b) for the NACA 64A010A and 
MBB-A3 airfoils, respectively. These results 
are compared with conventional method solutions 
to assess the accuracy of the tracking procedure 
when static twist is included. The results 
computed without static twist are those shown 
previously in Figs. 5(a) and 5(b). The overall 
effect of static twist on the flutter behavior 
of both airfoils Is an Increase in flutter speed 
index. Static aeroelastlc effects on flutter 
speed Index vary as the twist angles change with 



Fig. 7 Static aeroelastlc twist angle 

(“o - °r) f ° r ^e NACA 64A010A and 
MBB-A3 airfoils at flutter. 
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speed Index and Mach number. For example, at 
subcritical Mach numbers, static twist has 
little affect on flutter behavior, as would be 
expected based upon linear theory. As the Mach 
number Increases Into the transonic regime, the 
effects of static twist on flutter generally 
Increase. For the NACA 64A010A airfoil (Fig. 
8(a)), the flutter boundaries with and without 
static twist show only small differences between 
M = 0.70 and M = 0.74. At M « 0.76, the flow Is 
slightly supercritical and the static twist 
begins to have a larger effect on flutter. A 
shock forms on the airfoil upper surface between 
M = 0.76 and M = 0.78, which produces a large 
increase in flutter speed Index due to static 
twist. From M = 0.78 to M = 0.80, the static 


Conventional Tracking 

O Without static twist 

□ With static twist 



(a) NACA 64A010A airfoil with a? = 1.0°. 

Conventional Tracking 

O Without static twist 



(b) MBB-A3 airfoil with <v = -0.5°. 


Fig. 8 Comparison of flutter boundaries 

calculated with conventional and flutter 
boundary tracking procedures including 
static aeroelastlc twist. 


twist angle decreases and its effect on the 
flutter boundary is slightly reduced. For the 
MBB-A3 airfoil (Fig. 8(b)), there is a 

negligible difference In the flutter boundaries 
with and without static twist below M = 0.76, 
due to the subcritical nature of the flow. 
Between M = 0.76 and M = 0.78, the shock wave 
forms and static twist begins to affect the 
flutter boundary for the MBB-A3 airfoil in a 
manner similar to that for the NACA 64A010A 
airfoil. At M = 0.80, however, a more 
significant increase in flutter speed index 
occurs for the MBB-A3 airfoil due to the much 
larger (negative) static twist angle as shown in 
Fig. 7. A much larger value for flutter speed 
Index also results, since the shock wave on the 
MBB-A3 airfoil is much weaker in comparison with 
the NACA 64A010A airfoil. Edwards, et al. 7 
studied the effects of static twist for higher 
root angles of attack (a r = 2 to 5°) where 
they found more pronounced differences in the 
flutter boundaries of the two airfoils. The 
tracking procedure flutter boundaries including 
static twist in Figs. 8(a) and 8(b) are In 
excellent agreement with the respective 
conventional method solutions. The average 
errors for the boundaries with static twist are 
0.32% and 0.19% for the NACA 64A010A and MBB-A3 
airfoils, respectively. The differences in the 
flutter boundary results with and without static 
twist demonstrate the importance of static 
deformation on flutter behavior. The inclusion 
of static twist is therefore a requirement for 
performing realistic flutter analyses for 
wings. The tracking procedure thus contains all 
of the essential features for application to 
practical three-dimensional cases. 

Flutter Boundary Tracking Results as a Function 
o f Elastic Axis Location 

Flutter boundary tracking results as a 
function of elastic axis location are shown in 
Fig. 9 for the NACA 64A010A airfoil at H = 

0.80. Calculations were performed both with and 
without static aeroelastlc twist. The airfoil 
mass properties were held fixed by requiring 
that x a = -(a + 0.2) which fixes the mass 

center and that r 2 = r 2 + x 2 (where r = 0.49) 
a eg a eg 

which fixes the radius of gyration about the 
center of gravity (r C g). The step in elastic 
axis location was selected as &a = 0.2, and 
flutter boundaries were calculated for the range 
-2.0 ^ a < -0.6. The results for a = -2.0 are 
the same as those shown in Fig. 8(a) at M = 
0.80. As shown in Fig. 9, the flutter speed 
Index decreases as the elastic axis location is 
moved aft from a = -2.0 to a = -0.6. The 

decrease in speed index is due to a lower pitch 

mode frequency which increases the coupling 
between the plunge and pitch modes. For 

example, at a = -2.0, the wind-off coupled 
plunge and pitch natural frequencies are 71.33 
rad/sec and 535.65 rad/sec, respectively. At 
a = -0.6, these frequencies become 78.23 rad/sec 
and 165.26 rad/sec, respectively. The flutter 
boundary with static twist has flutter speed 
indices which are greater than those for the 
boundary without static twist. This increase in 
flutter speed index due to static twist is a 
function of the twist angle (ag - a,-). The 
twist angle decreases from (oq - “r) = -0.54° 
at a = -2.0 to (a,, - a,.) = -0.36° at 
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Conventional Tracking 

O Without static twist 

□ with static twist 



Fig. 9 Flutter boundaries as a function of 
elastic axis location for the NACA 
64A010A airfoil at M = 0.80. 


a = -0.6. Consequently, the Increase In flutter 
speed Index due to static twist Is greatest at 
a = -2.0. 


Concluding Remarks 

An automated flutter boundary tracking 
procedure has been developed for the efficient 
calculation of transonic flutter boundaries. 
This new procedure uses aeroelastlc responses 
computed with the XTRAN2L unsteady transonic 
small -disturbance code, to march along the 
boundary by taking steps In speed and Mach 
number. The flutter boundary tracking procedure 
therefore reduces the number of response 
calculations previously required to determine a 
transonic flutter boundary, and provides a 
complete boundary In a single job submission. 
Furthermore, the tracking procedure reduces 
computational costs since only two response 
calculations are required per Mach number. 

Flutter boundary results were presented for 
a simple aeroelastlc system consisting of a 
typical airfoil section oscillating with pitch 
and plunge degrees of freedom, to demonstrate 
the tracking procedure. These flutter 
boundaries were In good agreement with "exact" 
boundaries calculated using the conventional 
method. With the flutter boundary tracking 
procedure, the elapsed analysis time has been 
reduced from two to four weeks to one day 
turnaround, and the computational cost 
approximately halved. 

To be able to conceptually apply the 
tracking procedure to realistic 

three-dimensional cases, the procedure was 
extended to Include static aeroelastlc twist as 
a two-dimensional simulation of the static 
deformation of a wing. The tracking procedure 
flutter boundaries computed with static twist 


were in excellent agreement with the "exact" 
solution, and cost no more than the boundaries 
obtained without static twist. Therefore, the 
flutter boundary tracking procedure Is accurate 
and contains all of the essential features that 
are required to apply it to practical 
three-dimensional cases. 

Additionally, the tracking procedure is not 
constrained to determine flutter boundaries as a 
function of Mach number. The procedure is 
applicable to other parameters influencing the 
flutter boundary, making it potentially useful 
as a design tool. A sample calculation was 
presented showing a flutter boundary as a 
function of elastic axis location, thereby 
demonstrating this capability. 

Finally, the present work is considered to 
be an “initial" capability since flutter 
boundary tracking is a first attempt to provide 
an efficient automated procedure. Future work 
will be aimed at extensions of the present 
capability as well as other approaches of 
transonic flutter boundary determination. 
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